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Summary 


Following  the  previous  report,  where  a  literature  study  on  current  sensor  capacities  has  been 
carried  out  together  with  the  mathematical  definition  of  automatic  domain  splitting  (ADS)  and 
virtual  observatory  (VO),  this  work  contains  the  definition  of  the  initial  orbit  determination 
(IOD)  method  with  some  first  results,  a  literature  review  about  the  Manifold,  the  creation 
of  the  algorithms  and  architecture  for  the  ADS  and  the  initial  merging  of  the  ADS  into  the 
IOD  technique.  The  IOD  algorithm  is  explained,  implemented  and  comparisons  with  results 
from  literature  are  also  outlined.  In  specific  it  takes  as  input  optical  observations  and  gives 
as  output  the  state  of  the  object  at  the  central  time  of  observation  expanded  with  respect  to 
perturbations  on  the  observations.  Since  one  expansion  is  not  accurate  enough  to  describe  the 
full  initial  domain,  the  ADS  tool  is  introduced.  This  tool  is  constructed  on  the  mathematical 
definition  of  the  manifold.  The  initial  state  vector  is  then  defined  as  a  set  of  Charts  and  the 
propagator  integrates  in  time  this  initial  manifold  while  the  ADS  manages  the  convergence 
radius  of  manifold,  or  rather  the  convergence  radius  of  every  single  Chart.  Thus,  when  a  single 
truncated  power  series  (TPS)  is  not  sufficient  to  represent  the  whole  manifold,  that  is  when  the 
estimation  error  of  a  chart  is  bigger  than  the  threshold,  the  ADS  splits  the  manifold  and  goes 
ahead  with  the  propagation. 
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Chapter  1 

Introduction 


Every  day  thousands  of  detections  of  objects  orbiting  the  Earth  are  retrieved  by  observatories. 
However,  a  single  passage  above  the  station  is  not  sufficient  to  fully  determine  an  orbit.  The 
previous  report  showed  an  approach  to  perform  IOD  based  on  the  generation  of  an  AR  to  find  a 
space  of  solutions  rather  than  a  single  point.  Indeed,  due  to  the  observation  geometry  and  the 
uncertainty  related  to  sensor  accuracy,  timing  accuracy,  and  observer  state  knowledge,  it  is  not 
possible  to  obtain  a  single  orbital  state.  This  report  is  going  to  outline  a  new  method,  called  the 
Differential  Algebraic  Initial  Orbit  Determination  (DAIOD),  which  aims  to  describe  the  solution 
of  the  IOD  as  a  TPS  that  depends  on  variations  of  the  observations. 

Section  1.1  explains  the  output  of  a  typical  real  optical  observation.  Although  for  the  present 
work  the  observations  are  simulated  with  the  VO  (Section  2.2),  real  data  used  in  previous  works 
are  exploited  to  create  realistic  outputs.  Section  1.2  presents  the  relevant  nomenclature. 

1.1  The  observations 

A  typical  optical  observation  is  made  of  four  important  values: 

t  The  time  at  which  the  observation  is  made.  In  a  geostationary  Earth  orbit  (GEO),  the  time 
between  two  subsequent  observation  is  typically  around  2—3  min; 

oc  The  longitude  of  the  observed  object  with  respect  to  the  observatory; 

S  The  latitude  of  the  observed  object  with  respect  to  the  observatory; 

<jp  The  precision  of  the  observation  (usually  in  arcseconds).  Typical  values  for  the  precision 
vary  from  1",  being  the  most  accurate  observation,  to  5— 10",  being  the  least  accurate. 

For  the  work  at  hand,  an  optical  observatory  with  a  very  large  field  of  view  (FOV)  has  been 
assumed.  As  outlined  by  Milani  et  al.  (2004),  the  closer  the  observations  are  in  time,  the  more 
difficult  it  is  to  perform  IOD.  For  this  reason,  given  a  full  passage  over  the  station,  the  first, 
middle  and  last  observations  are  considered  in  this  work,  which  results  in  about  8  min  separation. 
The  observations  used  to  perform  IOD  are  simulated  with  the  VO,  presented  in  Section  2.2. 

1.2  Nomenclature 

In  this  Section  a  quick  overview  of  the  nomenclature  used  for  the  work  at  hand  and  in  previous 
papers  regarding  IOD  of  space  debris  is  given. 
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Observation.  As  outlined  in  Section  1.1,  an  observation  is  made  of  a  time  t,  two  angles  oc 
and  8  and  the  associated  precision  a >. 

Very  short  arc  (VSA).  It  is  a  sequence  of  N  observations  where  an  object  is  found  to  move 
Milani  et  al.  (2004).  Although  introduced  with  such  a  name,  in  Milani  and  Gronchi  (2009) 
and  Worthy  III  and  Holzinger  the  name  has  been  changed  to  too  short  arc  (TSA),  while 
other  authors  simply  call  it  short  arc  (SA)  (Fujimoto  and  Scheeres,  2012).  However,  all 
the  papers  agree  on  its  definition:  due  to  the  short  interval  between  each  detection,  these 
observations  do  not  allow  for  the  definition  of  a  track,  but  still  contain  useful  information 
about  the  object. 

Observations  set.  It  is  the  set  of  three  observations  used  to  perform  IOD,  generally  the 
first,  middle  and  last  of  a  VSA. 

Attributable.  Milani  et  al.  (2004)  defined  an  attributable  the  useful  information  that  could 
be  extracted  from  a  TSA.  It  is  made  of  two  angles  and  two  angular  rates: 

A  =  (oc,8,6c,8) 

where  generally  the  two  angles  coincide  with  the  middle  observation  and  the  angular 
rates  are  calculated  with  the  remaining  data  from  the  observation  set.  This  definition  is 
now  used  by  all  authors  to  describe  the  4  D  vector  containing  partial  information  about 
the  object  observed. 

(^,^)-domain.  Considering  the  variations  of  the  observed  angles  as  gaussians  with  zero 
mean  and  Op  •  variance 


f  5a.i  ~  N  (o,  oj^j 
jy,  ~v(o  ,cr^.) 


for  i  =  1,2,3 


(1.1) 


the  (a,  8 ) -domain  is  the  6 D  region  containing  the  3 crP  variation  of  the  angles. 

Orbit  set  (OS).  Given  the  state  of  the  object  as  a  function  of  8a  and  58,  the  orbit  set  (OS) 
is  the  set  of  the  orbital  states  enclosed  by  the  range  the  state  function  over  a  given  (a,  8)- 
domain.  Once  the  OS  is  defined,  it  is  possible  to  retrieve  other  values  as  a  function  the  state, 
thus  also  having  them  expanded  with  respect  to  the  same  (a,  8) -domain.  Functions  that 
will  be  used  later  are,  for  example,  the  range  and  the  range-rate:  p{5a,  58)  and  p(8a,  58). 

Admissible  region  (AR).  It  is  the  2D  plane  generated  by  the  two  degrees  of  freedom  of 
the  Attributable,  that  is  the  (p,p)  plane,  where  p  is  the  range  and  p  is  the  range-rate.  It  is 
the  region  where  the  attributable  places  the  information,  respecting  physical  constraints 
such  as  the  energy  law  and  minimum /maximum  distance  from  the  Earth.  Alternatively, 
it  can  be  described  as  the  range  of  the  p(5oc,  88)  and  p  ( doc,  58)  functions,  given  the  initial 
(a,  8) -domain. 
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Chapter  2 

SPICE  and  the  Virtual  observatory  (VO) 


This  chapter  introduces  two  important  tools  used  in  the  work  at  hand:  Spacecraft  Planet 
Instrument  C-matrix  Events  (SPICE)  and  the  VO.  The  first  one  is  a  powerful  tool  developed 
by  the  Navigation  and  Ancillary  Information  Facility  (NAIF)  that  is  able  to  assist  scientists  in 
planning  and  interpreting  scientific  observations  from  space-borne  instruments,  and  to  assist 
engineers  involved  in  modeling,  planning  and  executing  activities  needed  to  conduct  planetary 
exploration  missions.  The  VO  is  the  tool  developed  to  simulate  the  observations  and  is  based  on 
SPICE  functions. 

2.1  SPICE 

The  Navigation  and  Ancillary  Information  Facility  (NAIF),  acting  under  the  directions  of  Na¬ 
tional  Aeronautics  and  Space  Administration  (NASA)'s  Planetary  Science  Division,  has  built 
an  information  system  named  SPICE,  which  use  extends  from  mission  concept  development 
through  the  post-mission  data  analysis  phase,  including  help  with  correlation  of  individual 
instrument  data  sets  with  those  from  other  instruments  on  the  same  or  on  other  spacecraft 
(NAIF,  2016). 

As  the  name  reveals,  SPICE  contains  a  big  variety  of  information: 

S  Spacecraft  ephemeris.  When  this  file  is  available,  the  satellite  state  can  be  retrieved  at  any 
time  within  the  interval  of  definition  of  the  file. 

P  Planet,  satellite,  comet,  or  asteroid  ephemerides,  or  more  generally,  location  of  any  target 
body,  given  as  a  function  of  time.  It  also  includes  physical,  dynamical  and  cartographic 
constants  for  target  bodies,  such  as  size  and  shape  specifications,  and  orientation  of  the 
spin  axis  and  prime  meridian. 

I  Instrument  description  kernel,  containing  descriptive  data  peculiar  to  a  particular  scientific 
instrument,  such  as  field-of-view  size,  shape  and  orientation  parameters. 

C  Pointing  kernel,  containing  a  transformation,  traditionally  called  the  "C-matrix",  which 
provides  time-tagged  pointing  (orientation)  angles  for  a  spacecraft  bus  or  a  spacecraft 
structure  upon  which  science  instruments  are  mounted. 

E  Events  kernel,  summarizing  mission  activities  -  both  planned  and  unanticipated. 

The  Events  kernel  is  rarely  used.  For  the  DAIOD  algorithm,  only  the  S  and  P  kernels  will  be 
used.  In  particular,  they  will  be  necessary  to  create  the  observation  (as  explained  in  Section 
2.2),  retrieve  the  observatory  and  Earth  states  at  certain  times,  easily  recover  body  constants 
and  make  conversions  between  units  of  measurements.  While  the  Earth  ephemerides,  the  body 
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constants  and  the  conversion  files  are  already  available  in  the  SPICE  download  package,  the 
ephemerides  of  the  observatory  have  to  be  built  with  some  SPICE  toolkits:  after  loading  the 
Earth  ephemerides  file  and  the  file  containing  the  conversion  between  the  Earth  inertial  and 
Earth-fixed  reference  frames,  the  ephemerides  of  the  observatory  on  the  Earth  can  be  computed 
with  the  PINPOINT  and  MKSPK  toolkits. 

2.2  The  Virtual  observatory 

The  VO,  referred  to  as  IRIS  from  now  on,  is  a  tool  that  can  simulate  optical  and  radar  survey 
scenarios.  It  is  based  on  a  high  fidelity  numerical  propagator  (including  Earth  geopotential,  third 
body  perturbations,  solar  radiation  pressure  and  atmospheric  drag)  capable  of  accurately  predict¬ 
ing  the  motion  of  object  in  different  orbital  regimes:  low-Earth  orbit  (LEO),GEO,medium-Earth 
orbit  (MEO)  and  high-Earth  orbit  (HEO).  The  object  ephemerides  are  automatically  converted 
into  SPICE  kernels  (NAIF,  2016)  to  accurately  simulate  observations  from  any  ground  location. 
The  observations  are  simulated  implementing  different  visibility  constraints  and  producing 
outputs  compatible  with  real  sensors. 

The  IRIS  tool  is  based  on  the  automatic  generation  of  SPICE  kernels  for  both  space  object 
trajectories  and  observatories.  The  inputs  set  by  the  user  are  (see  yellow  boxes  in  Figure  2.1): 

•  object  initial  conditions; 

•  object  parameters; 

•  simulation  window; 

•  observer  location; 

•  instrument  type  and  accuracy  statistics. 

The  initial  conditions  of  the  space  object  can  be  given  in  different  formats,  i.e.,  two-line  elements 
(TLE),  Cartesian  coordinates,  or  osculating  classical  orbital  elements.  Once  the  user  has  set 
the  initial  space  object  state  and  the  model  parameters,  the  tool  computes  the  object  trajectory 
through  a  high-fidelity  orbital  propagator,  referred  to  as  Accurate  Integrator  for  Debris  Analysis 
(AIDA),  and  generates  its  corresponding  kernel  through  SPICE  functions.  Figure  2.1  shows  the 
IRIS  tool  together  with  SPICE  functionalities  and  AIDA. 


Figure  2.1:  IRIS  software  tool  description. 
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Once  the  object  and  observatories  kernels  are  built  with  the  necessary  toolkits,  all  the  power¬ 
ful  features  of  SPICE  can  be  exploited  to  study  the  relative  geometry  between  the  object  and 
observer.  This  includes  the  possibility  of  generating  observables  in  a  large  number  of  inertial,  lo¬ 
cal,  and  user  defined  reference  frames,  avoiding  all  the  nuisances  associated  to  time  conversions, 
reference  frame  conversion,  and  observation  details  such  light  time  and  stellar  aberration.  The 
input  deriving  from  the  sensor  type  is  used  to  quickly  compute  the  visibility  windows  and  to 
simulate  the  associated  observables.  More  specifically,  all  constraints  deriving  from  the  sensor 
type  are  evaluated  by  exploiting  the  functionalities  offered  by  SPICE.  For  instance,  whenever 
one  considers  an  optical  sensor,  some  aspects  such  as  the  sky  background  luminosity,  the  object 
illumination,  and  the  object  elevation  can  be  taken  into  account  to  define  the  observability 
windows.  Finally,  the  relative  geometry  between  the  observatories  and  observable  objects  can 
be  used  in  combination  with  user-defined  measurement  noises  or  as  input  for  sensor  simulators 
to  generate  simulated  observations. 

The  output  of  the  IRIS  tool  is  the  observations  set  defined  in  Section  1.2. 
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Chapter  3 

Initial  orbit  determination  (IOD) 


Orbit  determination  (OD)  refers  to  the  use  of  a  set  of  techniques  for  estimating  the  orbits  of 
objects  and  is  typically  divided  into  two  phases.  When  the  number  of  observations  is  equal  to 
the  number  of  unknowns,  a  nonlinear  system  of  equations  need  to  be  solved.  This  problem  is 
known  as  initial  orbit  determination  (IOD).  When  many  more  observations  are  taken  over  an 
orbit  arc  of  adequate  length,  accurate  orbit  determination  (AOD)  can  be  performed  (Armellin 
et  al.,  2015).  As  outlined  in  the  previous  report,  there  are  currently  two  different  types  of  sensors 
able  to  survey  the  sky:  radar  and  optical.  So  far,  the  latter  has  been  taken  into  account  to  build 
the  algorithm  for  IOD  and  the  type  of  observation  was  explained  in  Section  1.1.  Then,  Section 

2.1  outlined  some  SPICE  functions  necessary  to  retrieve  constants  and  fixed  objects  states.  Now, 
other  tools  necessary  to  complete  the  IOD  algorithm  are  explaine:  in  Section  3.1  some  differential 
algebra  (DA)  tools  are  introduced.  Section  3.2  shows  two  well-known  algorithms  and,  lastly, 
the  mathematical  definition  of  the  IOD  algorithm  will  be  performed  in  Section  3.3.  Results  and 
comparisons  with  literature  are  then  outlined  in  Section  5.1. 

3.1  Introduction  to  differential  algebra  (DA)  tools 

DA  supplies  the  tools  to  compute  the  derivatives  of  functions  within  a  computer  environment. 
More  specifically,  by  substituting  the  classical  implementation  of  real  algebra  with  the  imple¬ 
mentation  of  a  new  algebra  of  Taylor  polynomials,  any  function/  of  v  variables  is  expanded  into 
its  Taylor  polynomial  up  to  an  arbitrary  order  n  with  limited  computational  effort.  In  addition  to 
basic  algebraic  operations,  operations  for  differentiation  and  integration  can  be  easily  introduced 
in  the  algebra,  thus  finalizing  the  definition  of  the  differential  algebra  structure  of  DA  (Berz, 
1986, 1987).  In  the  following  subsections,  the  expansion  of  the  solution  of  parametric  implicit 
equations  and  the  nonlinear  mapping  of  the  estimate  statistics  are  explained.  The  descriptions 
are  taken  from  Armellin  et  al.  (2012)  and  Armellin  et  al.  (2015). 

3.1.1  Expansion  of  the  solution  of  parametric  implicit  equations 

Well-established  numerical  techniques  (e.g.,  Newton's  method)  exist,  which  can  effectively 
identify  the  solution  of  a  classical  implicit  equation 

/(*)  =  0.  (3.1) 

Suppose  an  explicit  dependence  on  a  parameter  p  can  be  highlighted  in  the  previous  function  /, 
which  leads  to  the  parametric  implicit  equation 

f(x,p)=  0.  (3.2) 
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Suppose  the  previous  equation  is  to  be  solved,  whose  solution  is  represented  by  the  function  x(p) 
returning  the  value  of  x  solving  (3.2)  for  any  value  of  the  parameter  p.  Thus,  the  dependence 
of  the  solution  of  the  implicit  equation  on  the  parameter  p  is  of  interest.  DA  techniques 
can  effectively  handle  the  previous  problem  by  identifying  the  function  x(p)  in  terms  of  its 
Taylor  expansion  with  respect  to  the  parameter  p.  The  DA-based  algorithm  is  presented  in 
the  following  for  the  solution  of  the  scalar  parametric  implicit  Eq.  (3.2);  the  generalization  to  a 
system  of  parametric  implicit  equations  is  straightforward. 

The  solution  of  (3.2)  is  sought,  where  sufficient  regularity  is  assumed  to  characterize  the 
function  /;  i.e.,  /  e  Ck+1.  This  means  that  x(p)  satisfying 

f(*(p),p)=  0  (3-3) 


is  to  be  identified.  The  first  step  is  to  consider  a  reference  value  of  the  parameter  p  and  to 
compute  the  solution  x  by  means  of  a  classical  numerical  method;  e.g.,  Newton's  method.  The 
variable  x  and  the  parameter  p  are  then  initialized  as  k- th  order  DA  variables,  i.e., 

fxl  =  x  +  Sx 

[Pi  =  p  +  Sp. 


A  DA-based  evaluation  of  the  function  /  in  (3.2)  delivers  the  k-th  order  expansion  of  /  with 
respect  to  x  and  p: 

Sf  =  Mf(Sx,5p),  (3.5) 

where  M  f  denotes  the  Taylor  map  for  /.  Note  that  the  map  (3.5)  is  origin-preserving  as  x  is  the 
solution  of  the  implicit  equation  for  the  nominal  value  of  the  parameter  p;  thus,  Sf  represents 
the  deviation  of  /  from  its  reference  value.  The  map  (3.5)  is  then  augmented  by  introducing  the 
map  corresponding  to  the  identity  function  on  p  (i.e.,  Sp  =  Tp(Sp))  ending  up  with 


'  Sf 

'  Mf  ' 

Sx 

Sp 

.  xv  . 

Sp 

The  k- th  order  map  (3.6)  is  inverted,  obtaining 


Sx 

'  Mf  ' 

-1 

'  Sf 

Sp 

XP 

Sp 

(3.6) 


(3.7) 


The  inversion  is  a  built-in  tool  in  the  C++  implementation  of  DA,  called 
compute  the  k- th  order  Taylor  expansion  of  the  solution  manifold  x(p) 
evaluated  for  Sf  =  0 


Sx 

'  Mf  ' 

-1 

0 

Sp 

xp 

Sp 

DACE.  As  the  goal  is  to 
of  (3.2),  the  map  (3.7)  is 


(3.8) 


The  first  row  of  map  (3.8) 


Sx  =  Mx(Sp), 


(3.9) 


expresses  how  a  variation  of  the  parameter  Sp  affects  the  solution  of  the  implicit  equation  as  a 
k- th  order  Taylor  polynomial.  In  particular,  by  plugging  map  (3.8)  in  first  of  (3.4)  we  obtain 


[x]  —  x  +  Mx(Sp),  (3.10) 

which  is  the  k- th  order  Taylor  expansion  of  the  solution  of  the  implicit  equation.  For  every  value 
of  Sp,  the  approximate  solution  of  fix,  p)  =  0  can  be  easily  computed  by  evaluating  the  Taylor 
polynomial  (3.10).  The  solution  obtained  by  means  of  map  (3.10)  is  a  Taylor  approximation  of 
the  exact  solution  of  Eq.  (3.2).  The  accuracy  of  the  approximation  depends  on  both  the  order  of 
the  Taylor  expansion  and  the  displacement  Sp  from  the  reference  value  of  the  parameter.  Thus, 
a  careful  analysis  is  always  mandatory  to  time  the  expansion  order  to  assure  that  Eq.  (3.10)  is 
sufficiently  accurate  for  the  entire  range  of  p  we  are  interested  in. 
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3.1.2  Nonlinear  mapping  of  the  estimate  statistics 

Consider  a  random  variable  x  e  3?”  with  probability  density  function  (PDF)  p(x)  and  a  second 
random  variable  y  e  1iv"f  related  to  x  through  the  nonlinear  transformation 

y  =  /(*)•  (3di) 

The  problem  is  to  calculate  a  consistent  estimate  of  the  main  cumulants  of  the  transformed  pdf 

p(y)- 


The  Taylor  expansion  of  y  with  respect  to  deviations  Sx  can  be  obtained  automatically  by 
initializing  the  independent  variable  as  a  DA  variable  and  evaluating  (3.11)  in  DA  framework. 
For  the  Z-th  component  of  y ,  this  procedure  delivers 

[y*]  =  /KM)  =  y«  +  Myi{Sx)  =  E  h,Vx...Vn  •  •  •  •  <$*»",  (3.12) 

Pi  4 - i“Pn^ 

where  in  this  expression  yz  is  the  zero-th  order  term  of  the  expansion  map,  and  are  the 

Taylor  coefficients  of  the  resulting  Taylor  polynomial 


The  evaluation  of  (3.12)  for  a  selected  value  of  Sx  supplies  the  Zc-th  order  Taylor  approximation 
of  \ji  corresponding  to  the  displaced  independent  variable.  The  Taylor  series  in  the  form  (3.12) 
can  be  used  to  efficiently  compute  the  propagated  statistics  (Valli  et  al.,  2012;  Park  and  Scheeres, 
2006).  The  method  consists  of  analytically  describing  the  statistics  of  the  solution  by  computing 
the  Z-th  moment  of  the  transformed  pdf  using  a  proper  form  of  the  Z-th  power  of  the  solution 
map  (3.12).  The  result  for  the  first  two  moments  is 


y*  =  £{  [Vi] }  =  E  ci/Pl...PnE{K  ■  •  •  Sx? } 

Pi  4 - \~Pn^k 

^ViVi  —  £{([y/]  “  m)(\Vj\  —  Vj))  —  2  Ci,Pl"-PnCj,qi'~qnE{^Xl 

Pi  4 - \-pn^k, 

Cjl~\ - h  Cjn^zk 


■■5xpnn+cln},  (3.14) 


where  are  the  Taylor  coefficients  of  the  Taylor  polynomial  describing  the  Z-th  component 

of  [y]  (in  the  covariance  matrix  formula,  the  coefficients  Cj/Pl^^n  and  Cj/Cjl^^n  are  updated  to 
include  the  subtraction  of  the  mean).  Note  that  the  expectation  values  on  the  right  side  of  Eq. 
(3.14)  are  function  of  the  known  p( x). 


When  x  is  a  Gaussian  random  variable,  its  statistics  are  completely  described  by  the  first 
two  moments,  i.e.  the  mean  y  and  the  covariance  matrix  P.  The  expectation  value  terms  of 
Eq.  (3.14)  are  thus  functions  of  the  initial  mean  and  covariance  only  and  they  can  be  computed 
applying  Isserlis's  formula  (Isserlis,  1918).  The  resulting  moments  are  then  used  to  describe  the 
transformed  PDF. 


3.2  Gauss's  and  Lambert's  Algorithms 

3.2.1  Gauss's  Algorithm 

Gauss's  algorithm  works  in  double  precision  and  takes  as  input  the  times  of  observations 
(hr  hr  h)r  the  positions  of  the  observatory  at  these  times  (R],  R2,  R3)  and  the  direction  cosine 
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vectors  (/91/92/&3).  The  algorithm  estimates  the  slant  ranges  (pi,  p2,  P3)  in  order  to  obtain  the 
object  positions 

rz  —  Rj  +  pi  p.,  where  i  —  1, 2, 3  (3.15) 

in  two-body  dynamics.  Figure  3.1  shows  the  geometry  of  the  problem.  The  description  of  this 
algorithm  refers  to  Curtis  (2015). 


Figure  3.1:  Geometry  of  input  and  output  values  for  Gauss's  algorithm.  Modified  from  Curtis 
(2015). 


Vector  equation  3.15  is  made  of  9  scalar  equations  and  12  unknowns.  Three  additional  scalar 
equations  can  be  added  by  remembering  that  the  three  position  vectors  must  lie  in  the  same 
plane,  due  to  the  angular  momentum  conservation.  Thus: 


r2  =  C\T\  +  C3T3.  (3.16) 

However  Equation  3.16  introduces  two  new  unknowns,  bringing  the  total  of  equations  and 
unknowns  to  respectively  12  and  14.  Another  consequence  of  the  two-body  dynamics  is  that  r 
and  v  can  be  expressed  in  terms  of  the  state  vector  at  any  given  time  by  means  of  the  Lagrange 
coefficients.  In  this  case  the  two  extremes  r\  and  7*3  can  be  defined  in  terms  of  the  state  at  central 
time  t2: 

|ri=/!r2  +  gl02  (317) 

(/3  =  hr2+g?,V2 

where  /„  g,  depend  on  f2,  i'2  with  reasonable  precision.  This  means  that  6  equations  and  4 
unknowns  are  added,  having  at  the  end  18  equations  and  18  unknowns,  to  be  solved  to  obtain 
the  values  0,7*2  and  r 3. 

The  algorithm  starts  by  calculating  quantities  that  are  going  to  be  useful  in  the  process: 


Vi  =  Pj  X  Pkr 

where  {i,j,k}  =  {1,2,3},  {2,3,1},  {3,1,2} 

(3.18) 

Tl  =  h  ~  h, 

t-H 

1 

CO 

II 

CN 

1 

CO 

H^k 

II 

(3.19) 

Djj  =  Ri  •  Pj, 

where  i  =  1, 2, 3 

(3.20) 

Do  =  h  ■  Pi 

(3.21) 
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Once  these  quantities  are  defined,  one  can  find  A,  B,  and  E : 


*-£(■ 

—  D12—  +  D22  +  032  —  ^ 

T  T  / 

(3.22) 

B=6D„ 

D12  (t32  -  T2)  ^  +  D32  (r2  -  T2)  ^ 

(3.23) 

E  =  R2  •  p2 

(3.24) 

After  a  lengthy  calculation  well  shown  in  Curtis  (2015),  we  are  left  with  an  8^-degree  polynomial 

that  solves  for  the  position  vector  at  central  time: 

r\  +  a  r\  +  b  r\  +  c  =  0 

(3.25) 

where 

a  —  —  (A2  +  2AE  +  R.I) 

(3.26) 

b  =  -2}iB(A  +  E) 

(3.27) 

c  =  -ji2B2 

(3.28) 

In  order  to  avoid  numerical  issues,  the  polynomial  is  scaled  by  d  —  so  that  the  new 

coefficients  are  Due  to  the  nature  of  our  variable,  only  positive  real  solutions  are 

accepted.  The  polynomial  is  solved  in  the  computer  environment  by  calculating  the  eigenvalues 
of  its  relative  companion  matrix.  Since  the  polynomial  has  been  scaled. 


Y 2 —  A^'  *  d 


(3.29) 


where  i  depends  on  the  number  of  acceptable  solutions  and  A;  are  the  positive  real  eigenvalues. 
If  more  than  one  solution  is  found,  they  need  to  be  analyzed:  there  may  be  a  degenerate  solution, 
a  solution  out  of  the  observatory  range  of  view  or  a  single  observation  could  create  two  or  more 
possible  solutions  for  the  first  guess  of  the  IOD.  So  far,  around  1000  observations  have  been 
analyzed  with  Gauss's  algorithm  by  the  team  and  they  always  produced  one  solution.  However, 
if  more  than  one  acceptable  solution  is  found.  Equations  from  3.30  to  3.34  have  to  be  carried  out 
for  each  acceptable  so  that  the  three  position  vectors  are  found  for  each  solution.  This  would, 
at  the  end  of  the  algorithm,  give  as  output  i  solutions  for  a  single  observation. 

f'  = 1  -  h? ur}  Sl  =  Tl~  b  (|)  (M0) 

/3  =  1“iVT|  g3  =  T3~b(?J  (3-31) 


Now,  only  p\,  p2  p3  are  missing  to  compute  the  position  vectors: 

r\  (6Di0r  -  6D0ot3  +  6D20ti)  +  Doofiif  -  D20}iTf  -  D0ofiT2T3  +  D20^t2ti 


P 1  = 

.  uB 

P2  -  A  +  T3- 


(D0ji t2t3  -  Do^rf  +  6  D0r3rf ) 


P3  = 


r\  (6D12T  -  6D02T3  +  6D22T1)  +  Dq2^t|  -  D22/2x  ~  Dq2^T2T3  +  D22^T2Ti 

DqTi  [jiT1  -  }irl  +  6 tf) 


(3.32) 

(3.33) 

(3.34) 


Having  the  position  vectors  of  the  observatory  from  SPICE,  the  cosine  directions  from  the 
observations  and  the  slant  ranges  from  this  algorithm,  the  position  vectors  for  the  object  observed 
can  be  found. 
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3.2.2  Lambert's  Algorithm 

Lambert's  algorithm  takes  as  input  two  position  vectors,  the  At  between  them  and  the  gravi¬ 
tational  parameter  and  gives  as  output  the  velocity  vectors.  The  geometry  of  the  problem  is 
shown  in  Figure  3.4.  The  algorithm  can  work  with  both  double  precision  and  DA  variables. 


Figure  3.2:  Geometry  of  input  and  output  values  for  Lambert's  algorithm.  Modified  from  Curtis 
(2015). 


The  algorithm  starts  by  fixing  the  sign  of  A 9  depending  on  the  type  of  orbit:  prograde  or 
retrograde.  The  two  equations  that  uniquely  define  A 9  are: 


A  9  —  acos 


-r2' 


Then,  depending  on  the  si 


t\  r2 

[hz  =  (*T  x  r2)z 


(3.35) 


Table  3.1:  Definition  of  A 9  depending  on  type  of  orbit 


Prograde  Retrograde 

o  o 

A  V 

N  N 

A0  =  acos  ( —  A9  —  2n  —  acos  ( —  —^] 

\  n  r2  )  \  n  r2  ) 

A8  —  2n  —  acos  ( —  A 9  —  acos  ( —  —^] 

V  nr2  )  \nr2  ) 

the  four  possible  transfers.  In  the  picture  tm  stands  for  ///.  Recalling  the  description  made  in 
Subsection  3.2.1,  once  (r\,V\)  are  defined,  the  state  of  the  object  can  be  retrieved  at  any  time 
through  the  Lagrange  coefficients.  This  means  that  finding  the  Lagrange  coefficients  solves 
Lambert's  problem  (Curtis,  2015).  However,  it  is  not  possible  to  solve  the  problem  analytically. 
In  the  past,  several  methods  were  created,  as  outlined  by  Vallado  and  McClain  (2001),  here  the 
C++  implementation  by  Izzo  (2014)  has  been  used,  after  updating  it  to  be  able  to  accept  both 
double  precision  and  DA  variables. 
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Figure  3.3:  Geometry  of  the  four  possible  transfers  depending  on  type  of  orbit.  Modified  from 
Vallado  and  McClain  (2001). 


3.3  Mathematical  definition  of  the  initial  orbit  determination  (IOD) 
algorithm 


The  IOD  algorithm  takes  as  input  the  observation  of  an  object  as  defined  in  Section  1.1  and  gives 
as  output  the  TPS  of  the  object  state  at  the  central  time  of  the  observation.  To  do  so,  an  initial 
estimate  of  the  object  position  at  t\,  ti  and  1 3  is  obtained  through  Gauss's  algorithm  (Subsection 
3.2.1)  in  double  precision.  The  observatory  state  is  retrieved  with  SPICE  (Section  2.1)  and  the 
direction  cosines  are  found  though  the  observation  angles: 


Pi  = 


COS  Si  COS  Oii 
cos  Sj  sin  oii 
sin  Si 


(3.36) 


where  i  =  1,  2,  3  refer  to  the  observation  times. 


At  this  point  an  estimate  for  the  position  vectors  relative  to  the  observation  times  is  avail¬ 
able  in  double  precision.  Now,  the  velocities  have  to  be  computed.  Lambert's  algorithm  takes  as 
input  two  position  vectors  and  the  At  between  them  and  gives  as  output  the  velocity  vectors. 
This  means  that  by  computing  Lambert's  algorithm  twice  (from  fi  to  G  and  from  £2  to  h)  one 
should  be  able  to  retrieve  the  three  state  vectors.  However,  Gauss's  algorithm  does  not  ensure 
that  the  three  estimated  vectors  define  one  unique  orbit,  thus  it  does  not  ensure  that  the  two 
velocity  vectors  found  at  t2  (V2  and  v2  )  coincide,  as  can  be  seen  in  Figure  3.4.  To  fix  this  problem 
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Figure  3.4:  Output  of  Lambert  algorithm  taking  as  input  the  output  from  Gauss's  algorithm 


and  obtain  the  expansion  of  the  state,  Lambert's  algorithm  will  be  used  twice:  the  first  time 
to  modify  p  to  ensure  that  and  the  second  one  to  expand  the  solution  (7*2,  V2)  with 

respect  to  the  observation  angles  variations.  The  first  step  allows  us  to  improve  the  estimation 
of  the  ranges  made  in  Gauss's  algorithm,  while  the  last  step  is  very  important,  because  it  allows 
us  to  analyze  the  variations  in  the  state  vectors  due  to  variations  in  the  observations  just  by 
means  of  function  evaluations.  Indeed,  although  the  observations  are  our  knwon  values,  they 
not  free  of  errors:  sensor  accuracy,  timing  accuracy  and  observer  state  knowledge  all  influence 
the  observation. 

For  the  first  usage  of  Lambert's  algorithm,  the  values  p\,  p2,  p3  are  initialized  as  DA  variables. 
Equation  3.37  shows  the  mathematical  and  computer  adapted  definition: 

r 

Pi, DA  —  Pi, Gauss  +  dpi 
S  Pi, DA  —  Pi, Gauss  +  dpi 
P3,DA  =  P3, Gauss  +  dp3 

The  importance  of  counting  the  number  of  DA  variables  will  become  clear  later.  In  this  way,  the 
output  of  Lambert's  algorithm  are  the  velocity  functions  depending  on  variations  of  the  slant 
ranges.  In  particular: 

v~  =  z>"  (Spi,  Sp2)  (3.38) 

v+  =  v+  (Sp2,  Sp3 )  (3.39) 

With  the  goal  of  solving  the  discontinuity  in  t2,  the  Av  between  the  left  and  right  velocities  is 
calculated: 

Av  —  v£  -  =  v  (Spi,  Sp2,  5pz)  (3.40) 

By  forcing  Av  =  0  one  wants  to  find  the  dp  necessary  to  obtain  it.  Newton  method  for  DA,  well 
explained  in  Section  3.1.1,  is  used  here.  Indeed,  the  function  /  is  Av,  the  parameters  x  are  p,  and 
the  variations  Sp  are  Spi.  Thus,  from 

Av  (pi,  p2,  P3)  fyi,  $Pi,  dp 3)  =  0  (3.41) 

one  obtains: 

r 

Pl,Ll  =  Pi, Gauss  3“  Api 

<  Pi, Li  =  pi, Gauss  +  Ap2  such  that  Av  (p1/L1,  p2,Ll,  p3,Ll)  0,  0,  0)  =  0  (3.42) 

P3,L1  =  P3, Gauss  +  &p3 


{Pi, DA  —  Pi, Gauss  +  DA(1) 

Pi, DA  —  Pi, Gauss  +  DA(2)  (3.37) 

P3,DA  =  P3, Gauss  +  DA(3) 


19 

DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


At  the  end  of  this  step,  one  has  obtained  the  states  of  the  object  that  satisfy  the  constraint  of 
pertaining  to  a  unique  orbit.  However,  the  solution  is  expanded  with  respect  of  the  slant  ranges. 
This  is  not  useful,  since  one  wants  the  solution  in  terms  of  the  observations  and  not  in  terms  of 
the  output  itself.  For  this  reason,  Lambert's  algorithm  is  used  again.  This  time,  the  DA  variables 
initialized  are  the  observations.  Depending  on  the  total  number  of  DA  variables  allowed,  two 
methods  are  possible:  the  9-variables  and  6-variables  method.  The  first  one  can  exploit  again 
Newton's  method,  while  the  second  one  does  not  add  3  redundant  variables. 

•  9-variables  method. 

This  method  adds  6  DA  variables  to  the  three  already  existing,  which  are  re-initialized.  The 
non-constant  parts  of  the  angle  polynomials  are  scaled  by  the  precision  of  the  observation 
<Tp/(,  so  that  the  variation  of  the  final  solution  will  act  according  to  the  reliability  of  the 
observations.  Equation  3.43  shows  the  mathematical  and  computer  adapted  definition: 


Pi, DA  —  Pi, LI  +  $Pl 
Pl,DA  —  Pl,Ll  +  8p2 
P3,DA  =  p3,Ll  +  <$P3 
OCl ,DA  =  &l,obs  +  <7P,1  <5oii 
(X-2 ,DA  =  <X-2,obs  +  Vp,2  80C2 
0C3 ,DA  =  <X-3,obs  +  &P,3  8a3 
$1  ,DA  —  Si ,obs  +  aP,l  88\ 
$2, DA  —  $2,obs  +  aP,2  88 2 
,^3, ,DA  —  83 fibs  +  aP, 3  883 


Pi, DA  —  Pl,Ll  +  DA(1) 
P2,DA  —  P2,L1  +  DA(2) 
P3,DA  =  P3,L1  +  DA(3) 

OCl, DA  =  &l,obs  +  Pp,l  DA(4) 
OCl, DA  =  0i2,obs  +  Cp, 2  DA(5) 
0C3,DA  =  0C3,obs  +  °~P,3  DA(6) 
$1  ,DA  —  $l,obs  +  aP,l  DA(7) 
^2, DA  —  <$2,obs  +  aP,2  DA(8) 

„  3 3, DA  —  &3,obs  +  aP,3  DA(9) 


(3.43) 


Hence,  9  DA  variables  are  used.  This  means  that  at  the  beginning  of  the  algorithm  the 
three  states  have  the  following  dependencies: 


*1  =*i  (Spi,  5/x.i,  55\)  (3.44) 

Xl  —  *2  {Sp2>  $0i2,  S82)  (3.45) 

x3  =  *3  (Sp3,  80L3,  553)  (3.46) 

(3.47) 


After  applying  Lambert's  method  to  the  couples  [x\,  Xi)  and  (Xi,  x3),  the  constant  part 
of  Ay  is  null,  thanks  to  the  previous  modification  of  p  and  depends  on  all  9  variables. 
The  task  is  now  to  find  the  variation  in  the  slant  ranges  as  a  function  of  the  variation  in 
the  observations.  This  is  actually  fairly  easy,  since  Newton's  method  can  be  used  again. 
Indeed,  the  following  partial  polynomial  inversion  is  possible 

Av  ( Sp ;  Set,  86)=  0  =>  8p  =  8p  (8a,  83)  (3.48) 


thus  obtaining: 

r 

[Pl.Ll]  =  Pi, LI  +  Spi'DA  (80C,  83) 

"  [p2,L2]  —  p2,Li  +  8p2,DA  (8a,  83)  (3.49) 

^  [p3,L2\  —  P3,L1  +  8p3,DA  (8a,  83) 

Although  the  method  is  straightforward  and  fast,  at  the  end  of  the  algorithm  the  first  three 
DA  variables  are  left  unused.  This  is  not  a  problem  per  se,  but  it  wastes  memory  and  may 
create  instabilities  further  on.  For  this  reason,  the  6- variables  method  is  introduced. 
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•  6-variables  method. 

This  method  deletes  the  first  three  DA  variables  and  overwrites  them  with  the  6  DA  angles 
defined  in  the  9-variables  method.  Equation  3.50  shows  the  mathematical  and  computer 
adapted  definition: 


ft  1  ,DA  =  «l,ofcs  +  Pp,l 
<X-2, ,DA  =  <X-2,obs  +  &P, 2  2 

«3 ,DA  =  «3 ,obs  +  aP, 3  <$a3 
$1,DA  —  ,obs  +  VP,  1 
&2,DA  —  $2,obs  +  ^,2 
^3,DA  =  $3,obs  +  &P, 3 


«1,DA  =  «l,ofcs  +  <7p,l  DA(1) 
«2,DA  =  &2,obs  +  <7p,2  DA(2) 
a3,DA  =  «3,ofc  +  aP,3  DA(3) 
$1,DA  —  $l,obs  +  VP,  1  DA(4) 
&2,DA  —  ^2,obs  +  ^P,2  DA(5) 
^3, DA  —  h,obs  +  CrP/3DA(6) 


(3.50) 


At  this  point  p  is  the  output  of  the  first  Lambert's  algorithm  in  double  precision:  p  —  pi  \. 
The  mathematical  ground  of  this  method  is  the  first  order  Newton: 

Az>(p)=0  =>  pi+1=pi-J-^po)Av(pi)  (3.51) 

Here,  the  assumption  is  made  that  the  Jacobian  does  not  change  in  the  loop.  The  iteration  is 
carried  out  until  i  —  MaxOrder,  thus  until  the  highest  order  of  the  DA  variable  is  reached. 

Both  methods  mathematically  deliver  the  same  output:  the  slant  ranges  at  all  three  times  of  ob¬ 
servation  in  terms  of  the  variations  of  the  observations.  The  position  vectors  are  then  obtainable 
with  Equation  3.15,  while  the  velocities  can  be  calculated  with  one  last  Lambert's  procedure. 
The  only  difference  in  the  output  is  that  with  the  9-variable  method,  the  output  depends  on  9 
variables,  of  which  the  first  three  are  always  zero. 


An  important  outcome  of  this  method  is  that  one  not  only  obtains  the  point  solution,  but 
can  also  easily  calculate  variations  by  means  of  functions  evaluations.  The  original  domain  of 
the  variations  is  the  {pc,  S )  -domain,  while  the  range  of  this  domain  through  the  state  function 
(r2,v2)  is  called  the  orbit  set  (OS),  both  described  in  Section  1.2. 


The  solution,  however,  is  not  valid  for  any  variation  of  the  input.  Indeed  a  TPS  is  only  valid 
within  a  certain  region,  defined  by  the  radius  of  convergence.  Outside  this  convergence  disk,  the 
error  of  the  polynomial  approximation  is  bigger  than  a  certain  threshold,  which  means  that  the 
approximation  is  not  good  enough.  This  implies  that  one  single  TPS  is  not  able  to  describe  the 
entire  OS  and  the  initial  {pc,  5)- domain  needs  to  be  split  in  sub-regions.  With  this  division,  there 
would  be  as  many  TPSs  as  sub-regions,  keeping  the  error  of  the  approximation  below  a  certain 
threshold  for  the  whole  range  of  solutions.  The  tool  that  is  able  to  estimate  the  error,  divide  the 
domain  and  find  the  new  TPS  for  the  OS  is  called  the  ADS  tool  and  is  introduced  in  Chapter  4. 
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Chapter  4 

Automatic  Domain  Splitting 


Section  3.3,  ended  by  outlining  the  importance  of  analyzing  the  region  of  validity  of  a  TPS, 
since  it  defines  the  validity  of  a  DA  map.  To  solve  this  problem,  the  idea  is  to  create  a  tool 
able  to  calculate  the  limits  of  validity  of  a  TPS  and  able  to  split  the  initial  domain  in  two  or 
more  subsets,  when  a  single  Taylor  expansion  is  not  enough  to  represent  the  DA  map.  The  tool 
without  much  extra  computation  evaluates  a  new  Taylor  expansion  for  each  of  the  new  domains. 
This  approach  generates  a  list  of  domains  and  respective  TPSs,  whose  union  corresponds  to 
the  initial  DA  set.  In  this  way  all  the  TPSs  obtained  lie  on  a  region  within  their  convergence 
disk  and  the  approximation  error  can  be  controlled  and  kept  below  fixed  threshold.  This  can  be 
done  when  the  DA  map  is  created  as  well  as  when  it  is  propagated  in  time,  since  the  error  of  the 
Taylor  expansion  grows  along  the  trajectory  due  to,  for  example,  nonlinearities  in  the  dynamic. 
The  efficiency  of  DA  for  uncertainties  propagation  is  highlighted  in  Armellin  et  al.  (2010). 

The  mathematical  basis  on  which  the  ADS  operates  for  sub-domains  generation  and  structure  is 
described  by  the  concept  of  manifold  introduced  in  Section  4.1  and  4.2.  Then,  the  implementation 
of  the  ADS  routine  is  described  in  Section  4.3.  Sections  4.4, 4.5  and  4.6  describe  the  implemen¬ 
tation  of  the  ADS  algorithm  on  domains  belonging  to  different  Euclidean  spaces.  The  latter 
describes  the  space  where  the  IOD  algorithm  is  defined.  Lastly,  Section  4.7  introduces  the  ADS 
algorithm  for  propagation. 


4.1  Mathematical  Background  and  Definition 

Following  these  considerations  and  looking  for  a  new  advanced  method  to  represent  a  DA  set, 
the  mathematical  concept  of  Manifold  is  used  to  find  new  tools  for  the  work  at  hand.  To  better 
understand  the  concept,  some  mathematical  definitions  are  here  outlined  (Wittig,  2016) . 

Definition  1.  The  //-dimensional  manifold  Mn  is  a  topological  space.  Each  point  p  e  Mn  is 
contained  in  an  open  set  Up,  which  is  a  homeomorphism  of  an  open  set  of  the  Euclidean  space 
Rn. 

Using  this  definition  it  is  possible  to  define,  in  case  of  astrodynamics,  a  phase  space  as  a  6- 
dimensional  manifold.  More  detailed  properties  of  the  topological  manifold  are  presented  in 
many  available  Math  text  books,  such  as  Jeffrey  M.Lee  (2009). 

Definition  2.  A  Chart  is  the  couple  ( U,cp ),  where  U  c  M"  is  an  open  subset  of  the  manifold  and 


22 

DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


cp  :  U  1 — >  is  a  homeomorphism  of  U  into  Cartesian  space. 

Definition  3.  An  Atlas  si  is  a  collection  of  charts  (14,  <pu)  that  the  U a  form  a  open  cover  of  the 
manifold,  i.e.  [jaEA  ^  —  M. 

Every  topological  manifold  admits  such  an  atlas,  and  if  all  the  cpa  are  r-times  continuously 
differentiable,  we  refer  to  the  atlas  as  a  Cr-atlas.  The  atlas  is  not  unique,  since  many  different 
atlases  can  describe  the  same  manifold. 

To  summarize,  a  topological  n-dimensional  manifold  can  be  represented  by  an  atlas,  which  is 
composed  of  charts.  While  infinite  in  general,  many  manifolds  of  relevance  in  practice  (e.g. 
closed  manifolds)  can  be  described  by  a  finite  atlas. 

The  previous  definitions  are  similar  to  the  process  of  mapping  the  globe:  keeping  in  mind  the 
projection  of  a  part  of  the  globe  on  a  2-D  map,  the  part  of  the  globe  can  be  identified  with 
a  subset  of  the  manifold  and  once  it  is  defined  the  projection  function,  the  globe  section  can 
be  mapped  on  a  2-D  Cartesian  space.  Thus,  a  section  of  the  globe  and  its  projection  function 
constitute  a  Chart  for  the  globe.  Furthermore  does  not  exist  a  projection  function  able  to  map  the 
whole  globe  onto  a  single  one  Cartesian  space  chart  if  the  mapping  is  to  be  homeomorphic. 

4.2  Differential  algebra  (DA)  Manifold  Representation 

Since  the  q>a  functions  of  the  charts  are  bijective  functions,  it  is  also  possible  to  define  an  atlas 
using  the  inverse  maps,  therefore  the  inverse  charts  can  be  defined  as  (Va,  <pa_1)aeyl  where 

Applying  the  ideas  behind  to  the  DA  set  representation,  some  analogies  can  be  drawn.  A 
differential  algebraic  vector  (DAvector)  associated  with  a  fixed  domain  of  (—1, 1) v  can  be  seen 
as  representing  an  inverse  chart.  However,  just  like  it  is  impossible  to  map  the  whole  globe  with 
only  one  single  chart,  it  is  usually  not  possible  to  represent  a  DA  manifold  or  set  with  one  single 
DAvector.  Instead,  an  atlas  of  several  DA  vectors  is  used  to  cover  the  entire  DA  manifold 
The  choice  for  inverse  charts  mainly  comes  from  the  fact  that  it  is  much  easier  to  specify  a 
subset  of  Euclidean  space  for  the  inverse  chart  to  act  on  than  to  specify  a  subset  of  an  arbitrary 
manifold.  Practically,  this  means  that  the  DA  manifold  representation  of  in  our  case  the  TPS 
is  made  from  DAvector  mapping  from  Euclidean  space  to  phase  space.  Thus,  cpa~l  is  a  w- 
dimensional  DAvector  with  v  variables,  while  the  Euclidean  domain  for  simplicity  is  chosen  to 
be  (—1,1  )v.  That  definition  of  DA  Manifold  concides  with  that  given  by  Wittig  (2016). 

Definition  4.  A  p-dimensional  DA  manifold  M  embedded  in  a  zf-dimensional  space  is  defined 
by  its  finite  DA  atlas  si  of  DA  charts  ((—1, 1)4  (poc~l)  with  (p a_1  e  WDV. 

Moreover,  the  employment  of  inverse  maps  and  the  scaled  Cartesian  domain  lead  to  another 
simplification:  once  the  domain  is  fixed  to  (— 1,  l)v  it  does  not  need  to  be  stored  in  the  DA  chart 
to  save  memory. 

4.2.1  The  automatic  domain  splitting  (ADS)  and  differential  algebra  (DA)  mani¬ 
fold 

As  described  in  Chapter  3,  the  initial  state  is  generated  by  mapping  a  Cartesian  domain  and 
depends  on  the  nominal  values  and  their  uncertainties.  The  determination  of  a  Cartesian 
domain  and  its  maps  is  not  unique,  thus,  referring  to  the  manifold  analogy,  it  is  possible  to 
obtain  different  atlas  structures  which  represent  the  initial  state. 

For  the  DA  representation  of  the  DA  chart  and  the  DA  manifold,  Taylor  theory  is  used:  the  maps 
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are  Taylor  expansions  around  a  nominal  point,  that  is  a  map  representation  of  the  initial  state 
with  its  uncertainties.  This  mapped  region,  the  orbit  set  (OS),  can  be  propagated  to  a  chosen 
epoch.  For  the  Initial  state  determination,  the  map  is  a  function  that  goes  from  the  Cartesian 
space  to  phase  space,  while  for  the  state  propagation  the  map  goes  from  the  phase  space  to 
phase  space.  For  both  the  maps  defined  above,  the  convergence  radius  needs  to  be  respected. 
Indeed,  as  already  pointed  out,  the  estimation  of  the  error  between  the  true  value  and  the  map 
evaluation  is  strictly  linked  to  the  convergence  radius:  only  when  the  mapped  point  lies  within 
the  convergence  disk  it  well  approximates  the  real  state.  This  means  that  the  estimate  and  control 
of  the  error  is  a  crucial  point  in  our  work.  Since  the  ADS  provides  the  estimate  and  control  tools 
for  the  mapping  error,  it  is  briefly  introduced  here.Firstly,  the  maps  are  mathematically  defined. 
For  the  initial  state  determination  the  maps  goes  from  the  Cartesian  domain  U  e  Rv  to  the 
whole  DA  manifold,  which  represent  the  initial  state.  Thus,  the  atlas  representing  the  manifold 
and  defined  by  one  chart  is:  srf  —  {(IT, 7)}  where  7  :  (—1,1  )v  e  Rv  7((— 1,  l)v)  e  Rw.  The 
propagation  map  goes  from  the  phase  space  Rv  to  Rv  itself.  Indeed,  when  propagating  from 
Xo  to  xn,  the  function  associated  can  be  written  as:  xn  :  (—l/l)veRv^xn((—l/l)v)eRv.  For 
both  cases,  it  is  used  the  convergence  estimation  order  by  order  analysis  of  TPS  to  check  the 
estimation  error  of  the  map  (for  the  propagation  it  is  checked  at  every  time-step  tn)  and  if  this 
error  is  bigger  than  a  fixed  threshold,  the  ADS  tool  splits  the  initial  domain  U  into  two  or  more 
sub-domains  and  re-evaluates  the  maps  over  each  of  the  new  domains.  In  this  way  a  new  atlas 
that  represents  the  initial  manifold  is  obtained  and  can  be  written  as:  si  —  { (14, 7*)  }ocea  —  M 
where  |JaeA  14  =  U  and  7 dc:UoceRv  ^  7^(14)  e  Rw. 

4.3  Structure  and  Architecture  of  Algorithm 

The  mathematical  concepts  of  manifold  and  charts  need  an  implementation  for  computer  rep¬ 
resentation.  Three  main  classes  are  here  introduced  to  represent  each  component  of  the  DA 
manifold: 


Patch.  The  Patch  assumes  the  role  of  the  chart,  and  thus  the  scaled  Cartesian  domain  fixed 
with  (— 1,  l)v  and  the  map  onto  the  the  phase  space  are  stored  here.  To  optimally  store  the 
charts,  the  Taylor  expansion  of  the  map  and  SplittingHistory  are  stored  here,  since  the 
latter  implicitly  corresponds  to  the  domain.  Synthetically,  this  class  embeds  the  propriety 
of  a  DAvector  to  store  the  map  and  the  ID  SplittingHistory  to  identify  the  respective 
domain.  Lastly,  some  additional  member  functions  are  created  within  this  class  to  evaluate 
the  estimation  error  of  the  expansion  map,  to  establish  the  direction  of  the  split  and  to 
obtain  the  new  Patch  during  the  split. 

SplittingHistory.  As  outlined  before,  the  scaled  Cartesian  domain  and  inverse  maps  will  be 
stored  during  the  split.  However,  when  the  domain  is  split,  one  has  to  know  where  the 
new  domain  is  with  respect  to  the  initial  one.  For  this  reason  the  SpittingHistory  class  is 
introduced:  it  keeps  track  of  all  the  operations  performed  on  the  domain  through  a  vector 
that  stores  integer  values  representing  the  variable  direction  being  split.  To  complete 
the  class,  it  is  equipped  with  a  member  function  that  extrapolates  information  about  the 
domain  and  can  replay  it  after  or  during  the  splitting. 

Manifold.  The  Manifold  assumes  the  role  of  the  mathematical  manifold:  it  is  described  by  an 
atlas  that  is  the  union  of  several  Patches.  When  the  manifold  is  performed  by  the  ADS 
algorithm,  all  the  Patches  lie  within  the  convergence  disk,  thanks  to  the  construction  of  an 
error  controlling  function.  In  this  way,  the  ADS  algorithm  becomes  a  member  function 
to  decide  when  the  manifold  convergence  analysis  is  acceptable  or  not  acceptable  during 
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the  propagation.  In  short,  the  Manifold  is  defined  as  a  list  of  Patches  and  embeds  some 
member  functions  to  analyse  it,  to  extrapolate  additional  information  and  to  propagate  it. 

Hereafter  these  terms  will  be  used  during  the  explanation  of  the  work  and  to  give  results. 
Chapter  3  already  showed  some  useful  tools  provided  by  DA  to  solve  the  IOD  problem.  Starting 
from  these  results,  it  has  been  decided  to  build  a  new  algorithm  process  of  the  ADS,  based 
on  the  SplittingHistory,  Patch  and  Manifold  classes.  The  most  crucial  problem  with  the  DA 
representation  and  the  Taylor  expansion  is  the  calculation  of  the  error  between  the  real  function 
and  the  TPS  that  approximates  it,  thus,  the  theory  behind  the  error  estimation  function,  given  by 
Armellin  et  al.  (2010),  developed  within  the  ADS  algorithm  is  now  explained. 

Let  /  be  a  n  +  1  differentiable  function  /  e  Cn+1  and  let  its  Taylor  expansion  Pf  be  of  n  order, 
then  the  following  holds  for  the  error  introduced  by  the  expansion  around  the  origin: 

\f(Sx)  -Pf\  <  C-Sxn+1  (4.1) 

Now  let  er  be  the  maximum  error  of  the  expansion  Pf  on  the  domain  of  radius  r  >  0,  for  some  C 
that  depends  on  the  maximum  size  of  8  x.  If  we  reduce  the  domain  to  one  of  radius  \  also  the 
error  will  decrease  by  a  factor 

I/O**)  -Pf\<C  -Sxn+1  <  C  •  g)"+1  =  2^1  (4.2) 

This  justified  why  reducing  the  size  of  the  domain  improves  the  convergence  radius  of  the 
expansion  to  the  n+1  power. 

To  estimate  the  error  for  any  expansion  order,  the  coefficients  of  the  same  order  i  are  used  and 
the  sum  of  their  magnitude  S;  is  considered: 

Si  =  2  Kl  (4.3) 

I  *M 

When  variables  are  scaled  by  some  sufficiently  small  factor,  the  terms  of  any  convergent  power 
series  converge  at  least  exponentially.  When  this  is  the  case,  the  function  surely  converges  for  all 
arguments  smaller  than  1.  We  therefore  consider  this  the  convergence  radius  of  the  function  and 
the  exponential  decay  of  the  coefficients  by  order  is  the  feature  we  use  to  measure  convergence. 
An  exponential  fit  is  performed  is  performed  to  compute  the  parameters  A  and  B  ,  to  match  S; 
to  the  function 

Si  —  f(i)  —  A  •  exp{B  n )  (4.4) 

The  fit  function  /(z)  is  used  to  compute  the  value  f(n  +  1)  and  the  size  Sn+\.  If  the  Sn+\  is  too 
big  with  respect  to  the  error  threshold,  the  domain  is  split  and  new  expansions  are  evaluated 
over  the  two  new  domains,  repeating  this  process  sufficiently  often  ensures  that  the  convergence 
criterion  (i.e.  exponential  decay)  is  eventually  met. 

Once  the  error  estimation  is  established,  the  ADS  algorithm  can  work  on  different  domains  and 
maps  by  just  defining  the  initial  DAmanifold. 

The  common  procedure  for  different  types  of  map  functions  and  domains  is  now  described. 
The  first  step  is  the  definition  of  the  initial  domain,  thus  an  initial  Patch  and  lastly  the  Manifold 
through  the  DAmanifold  structure.  However,  this  first  Manifold  is  a  trial  manifold  because  it 
is  represented  by  a  mathematically  incorrect  atlas  structure.  Performing  the  ADS  algorithm 
which  automatically  analyses  the  manifold  and  then  the  ADS  constructs  the  atlas  structure  that 
represents  it.  An  overview  of  ADS  routine  and  the  connection  between  the  different  structures 
are  represented  in  Figure4.1  .  The  ADS  algorithm  needs  the  definition  of  a  function  that  maps 
the  Cartesian  space  onto  the  manifold.  This  function  can  be  defined  by  the  user  but  it  needs  to 
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Figure  4.1:  Overview  Framework  ADS. 


have  a  DAvector  structure,  which  means  that  the  input  and  output  must  be  a  DAvector,  for  the 
algorithm  to  be  able  to  store  it  in  a  DAPatch  and  then  into  a  DAManifold. 

This  ADS  algorithm  can  be  used  with  scalar  functions,  multi-variable  functions  and  ordinary 
differential  equations  (ODEs),  which  are  described  in  the  following  sections. 

4.4  Single- Variable  Scalar  Functions 

The  routine  of  the  ADS  algorithm  to  analyze  single-variable  scalar  functions  /  :  3ft  — >  3ft  is 
schematically  illustrated  in  Fig.4.2  on  page  28.  It  is  very  similar  to  the  process  previously 
described,  indeed  one  just  needs  to  define  the  initial  step  of  the  routine.  The  algorithm  starts 
from  the  initial  domain,  the  function  to  evaluate  the  scalar  value  takes  as  input  the  ID  Patch  of 
the  domain  and  returns  a  ID  Patch  as  output.  Once  the  Patch  of  scalar  function  is  obtained,  the 
ADS  algorithm  estimates  the  error:  if  the  error  is  below  a  fixed  threshold,  the  Patch  is  stored 
into  the  result  manifold,  otherwise  the  domain  is  split.  The  process  of  the  computation  of  the 
split  direction  gives  a  trivial  result,  that  is  the  single  variable  of  the  domain.  The  two  new  Patch 
sub-domains  are  inserted  into  the  trivial  manifold  and  the  process  is  resumed.  To  avoid  too 
small  subdomains  and  limit  the  runtime,  also  a  maximum  number  of  splits  Nmax  is  fixed.  This 
means  that,  it  can  happen  that  a  sub-domain  is  stored  in  the  results  Manifold  when  it  reaches 
the  minimum  allowed  size,  even  though  it  should  be  split  due  to  a  high  error  estimation. 

The  resulting  Manifold  contains  the  union  of  the  Patches,  which  represent  the  evaluation 
of  the  polynomials  in  each  of  these  sub-domains.  This  structure  provides  a  more  accurate 
approximation  than  the  one  achievable  with  a  single  polynomial  evaluation. 

4.5  Multi-variable  Scalar  Functions 

In  the  case  of  multi-variable  scalar  functions  /  :  — »  5ft,  the  ADS  routine  returns  a  ID  map 
and  always  analyzes  the  same  component.  After  the  computation  of  the  error  estimation  it 
evaluates  the  split  direction  as  shown  in  Figure  4.2.  In  this  case  the  result  is  not  trivial  because 
there  is  more  than  one  variable.  The  selection  of  the  optimal  split  direction  is  crucial  in  terms  of 
error  reduction.  This  direction  is  selected  through  the  procedure  used  in  the  propagation  of  the 
two-body  problem  by  Wittig  et  al.  (2015). 
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4.6  Multi-variable  Vector  Functions 


Also  for  the  case  of  multi-variable  vector  functions  /  :  ?ftk  — *  iv"' ,  the  routine  does  not  need 
any  adjustments  from  the  structure  represented  in  Fig.4.2  since  the  maps  are  Patch  elements. 
In  the  case  of  a  vector  valued  functions,  the  map  returned  by  the  routine  has  more  than  one 
component,  so  the  error  estimation  is  carried  out  for  each  dimension.  Then  the  component  of  the 
map  with  the  maximum  error  estimation  is  considered  and  the  polynomial  which  represents  this 
map  component  is  used  to  decide  the  splitting  direction.  The  only  caveat  in  the  vector  valued 
case  is  the  definition  of  the  meaning  of  "largest  error".  In  practice,  the  quantities  in  different 
components  of  the  vector  often  are  not  comparable  (e.g.  angles  and  distances),  and  hence  a 
single  splitting  threshold  is  not  appropriate.  In  such  cases,  a  weighting  factor  can  be  applied  to 
each  component  of  the  error  estimate  to  turn  them  into  comparable,  non-dimensional  units  of 
error.  If  the  maximum  residual  value  does  not  exceed  the  tolerance,  the  Patch  with  the  vector 
map  is  stored  in  the  result  manifold. 

The  multi-variable  vector  functions  are  very  useful  in  our  work  because  they  allow  us,  for 
example,  to  analyze  the  conversion  between  two  reference  frames  and  to  compute  the  evaluation 
of  the  initial  OS  determination  with  the  DAIOD  algorithm.  In  both  cases,  the  function  maps  two 
algebraic  structures  of  the  same  dimension  /  :  ^  Furthermore,  with  a  multi-variable 

vector  function,  one  can  also  map  the  AR  defined  in  Section  1.2. 


4.7  Flow  of  ordinary  differential  equations  (ODEs) 

The  routine  for  the  propagation  of  ODEs  is  described  in  Figure  4.3.  The  routine  starts  by 
considering  the  state  vector  at  the  initial  time  x(to)  —  Xo  and  it  is  stored  into  a  trial  manifold 
as  a  Patch  element.  In  the  propagation  case  the  Patch  needs  to  store  additional  information 
about  the  state  vector,  that  is  the  integration  step  which  identifies  the  propagation  time  as 
long  as  the  tolerance  and  state  vector  size  are  satisfied.  The  state  vector  is  propagated  until 
fnew  —  h  +  tstep  using  a  DA-based  integrator  to  obtain,  at  each  step,  the  Taylor  expansion  of  the 
solution  x  (neW/f )  with  respect  to  the  initial  conditions.  At  each  step,  the  tolerance  satisfaction  is 
checked.  If  the  tolerance  is  satisfied,  then  the  integrator  proceeds  to  the  next  time  step;  otherwise, 
the  integrator  stops,  the  domain  corresponding  to  the  current  state  X(0/z)  at  time  t\  is  split,  and 
two  new  state  vectors  are  generated  with  the  information  of  the  last  time  integration  t\.  They  are 
added  to  the  trial  manifold,  and  the  procedure  restarts  from  the  first  state  vector  in  the  list.  Every 
time  a  sub-domain  reaches  its  Nmax- th  split  or  the  final  simulation  time  tfinai,  the  set  is  saved  in 
the  result  manifold  with  the  corresponding  truncation  epoch.  The  integration  is  resumed  with 
a  new  state  vector  in  the  trial  manifold  until  the  list  is  empty,  which  means  that  the  routine  is 
completed. 
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Figure  4.2:  ADS  routine  to  sampling  the  scalar  and  vector  function  with  single  and  multi  variable 
domain 
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Figure  4.3:  ADS  routine  for  the  Taylor  expansion  of  the  differential  ordinary  equation 
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Chapter  5 

Results  and  Discussion 


5.1  IOD:  Admissible  region  (AR)  and  comparison  with  literature 

In  Section  1.2,  the  concepts  of  attributable  and  admissible  region  (AR)  are  introduced  as  part  of 
the  literature  available  on  IOD  of  space  debris.  Here,  they  are  represented  in  order  to  make 
comparisons  between  the  method  from  literature  (Milani  and  Gronchi,  2009)  and  the  current 
approach.  Figure  5.1  shows  the  AR  for  a  simulated  attributable.  The  energy  law  is  used  to 


Figure  5.1:  AR  from  literature.  The  green  area  is  the  CAR. 


constrain  the  region,  by  substituting  in  it  respectively  a  minimum/ maximum  semimajor  axis 
(a minr  Umax)  and  a  maximum  eccentricity  ( emax )•  In  this  way  three  constraining  inequalities  are 
obtained:  £amax  <  0,  £amin  <  0  and  £6max  <  0.  The  output  is  a  plot  with  an  implementation  of 
the  algorithm  by  Milani  and  Gronchi  (2009)  in  Matlab.  The  green  region  shows  all  the  possible 
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combinations  of  p  and  p  that  would  satisfy  the  physical  constraints  to  complete  the  attributable 
and  be  a  possible  state  for  the  observed  object.  As  can  be  seen,  it  is  a  fairly  big  region.  It  could  be 
shrunk  if  one  had  a  priori  information  about  the  object.  For  example,  knowing  that  the  satellite  is 
in  a  geostationary  transfer  orbit  (GTO)  would  put  more  stringent  constraints  on  the  eccentricity 
and  semimajor  axis,  thus  obtaining  the  CAR  in  Figure  5.2,  where  the  region  containing  circular 
orbits  has  been  canceled  out  and  only  eccentric  orbits  remain.  However,  a  priori  knowledge  is 


Figure  5.2:  AR  with  a  priori  knowledge.  The  green  area  is  the  CAR. 


in  general  not  available. 


The  AR  for  the  current  approach  is  defined  as  the  range  of  the  functions  p  and  p  starting 
from  the  (&,6) -domain.  Indeed,  being  the  equations  for  the  range  and  range-rate  a  function  of 
the  state,  p  and  p  also  depend  on  the  observations: 


p  =  Ikll  =  p(M) 

p  =  r-^-=p{u,5) 


(5.1) 


Subsection  3.1.2  described  how  to  map  the  PDF  of  a  variable  x,  which  in  this  case  are  See  and  88, 
into  a  transformed  PDF  through  a  non-linear  transformation,  that  is  Equation  5.1.  The  resulting 
box-shaped  AR  is  shown  in  Figure  5.3  superimposed  on  the  AR  presented  in  Figure  5.1.  For  this 
simulation,  high  precision  in  the  observation  was  assumed  (<Tp  —  1).  It  is  easy  to  check  that  the 
real  solution,  that  is  the  values  used  to  simulate  the  observation,  is  comprised  in  both  regions, 
but  the  box-shaped  region  obtained  with  the  approach  here  described  is  much  smaller:  indeed, 
it  is  already  possible  to  say  that  the  observed  object  is  in  a  HEO  with  a  small  or  zero  eccentricity. 
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dt/d,,  (km/s)  dt/d^  (km/s) 


Figure  5.3:  Comparison  between  ARs:  literature  vs.  presented  method.  Accurate  observation. 
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Figure  5.4  shows  the  comparison  of  ARs  for  a  less  accurate  observation,  while  Figure  5.5  shows 
the  comparison  for  an  accurate  observation  of  an  object  in  a  highly  elliptical  orbit. 


Figure  5.4:  Comparison  between  ARs:  literature  vs.  presented  method.  Inaccurate  observation. 


Figure  5.5:  Comparison  between  ARs:  literature  vs.  presented  method.  Highly  elliptical  orbit. 
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As  already  underlined  in  Section  3.3,  one  box  is  not  accurate  enough  to  describe  the  entire  region. 
Section  5.2.3. 2  will  show  the  same  ARs  presented  in  this  section  split  with  the  ADS  tool. 

Here,  a  quick  check  is  carried  out:  if  the  map  were  to  be  valid  everywhere,  the  error  between  the 
perturbed  observation  (referred  to  with  a  star  oc*)  and  the  angles  retrieved  with  an  evaluation  of 
the  perturbed  state  (referred  to  with  a  hat  &)  would  always  be  small  or  zero.  Figure  5.6  shows 
the  procedure 


Observation  Set  IOD 

(^5  ^5  ^5  ) 

+  5a 

+  SS  ? 

(«*.**) 


r  (5a,  55) 
v  (5a,  55) 

a  =  f(r,R) 
5  =  g(r,R) 

(a,  5) 


Figure  5.6:  Procedure  to  test  the  accuracy  of  the  map  (r,  v). 


where: 

*2=/(r,R)=atan2(^||,  jJ-|l)  (5.2) 

4  =  g(r,R)  =asin  (j|y-£|)  (5.3) 

r  is  the  object  position,  R  is  the  observatory  position  and  the  subscripts  refer  to  the  coordinate 
considered.  Figure  5.7  shows  the  contour  plot  of  the  procedure.  The  plots  represent  the  region 
of  the  (8a 2, 552)-plane  for  the  three  simulated  observations  already  used  before.  As  can  be  seen, 
the  map  is  valid  at  the  origin  (unperturbed  point),  but  a  deviation  of  few  arcseconds  would 
invalidate  it.  Remembering  that  observations  can  have  variations  <Tp  c  [1",  10"],  the  map  needs 
to  be  able  to  handle  them  without  introducing  significant  numerical  errors. 

5.2  Automatic  Domain  Splitting 

5.2.1  Single-Variable  Scalar  Function 

This  section  contains  the  first  step  to  validate  the  ADS  routine.  A  sine  function  has  been 
considered  to  check  the  accuracy  of  the  ADS  routine.  Using  the  guide  line  explained  in  Sect.  4.4, 
a  Patch  of  domain  x  e  [—3, 3]  has  been  created,  the  sine  function  has  been  defined  as  a  TPS  of 
order  5  and  the  tolerance  for  the  estimation  error  threshold  has  been  fixed  to  10-4.  Figure  5.8 
shows  the  result  for  the  approximation  without  (left)  and  with  (right)  the  ADS  routine. 

It  can  be  seen  that  8  sub-domains  are  sufficient  to  well  describe  the  function  by  respecting  the 
maximum  error. 

If  the  inverse-square  function  ^  is  used  as  map  function,  it  is  possible  to  appreciate  the  behavior 
of  ADS  routine.  In  Figure  5.9  the  inverse-square  function  is  plot  using  a  5-th  order  TPS,  10-4 
error  tolerance  and  fixing  the  domain  [—1.88,4.12]  (the  domain  should  be  such  that  it  does 
not  include  0  but  gets  close,  e.g.  [10-5, 3]  or  so  it  is  chosen  to  be  asymmetric  to  avoid  that  the 
mean  value  of  the  domain  or  any  of  the  subdomains  was  exactly  the  singular  point  which 
would  lead  to  a  failure  of  the  algorithm  as  the  Taylor  expansion  would  be  singular).  Lastly, 
the  maximum  number  of  split  for  each  sub-domain  is  fixed  to  7,  to  avoid  an  excessive  number 
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Accurate  GEO  observation 


Inaccurate  GEO  observation 
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-10  -5  0  5  10 
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Figure  5.7:  Logarithmic  contour  plot  to  check  the  validity  of  the  map  depending  on  perturbations 
of  &2  and  $2-  The  errors  are  shown  up  to  1",  then  colors  are  deleted. 


of  sub-domains.  The  left  side  of  Figure  5.9  shows  the  single  domain  expansion  vs.  the  exact 
solution,  demonstrating  that  a  single  expansion  is  not  enough  to  map  the  full  domain.  The 
right  side  shows  the  inverse-square  function  with  a  Manifold  of  34  Patches  which  increase  the 
precision  of  the  map.  The  same  figure  displays  on  the  x-axis  the  bound  of  the  subdomains: 
and  that  the  subdomains  of  minimum  size  are  close  to  the  singular  point.  The  reason  for  the 
presence  of  so  many  sub-domains  close  to  the  singular  point  is  that  the  Taylor  approximation  is 
poor  because  of  the  singularity,  thus  the  ADS  routine  keeps  splitting  the  domain  close  to  the 
singularity  point  until  the  maximum  number  of  splits  is  reached. 


5.2.2  Multi- Variable  Scalar  Function 

Carrying  on  the  validation  of  the  ADS  routine,  this  section  focuses  on  its  application  to  the 
multi-variable  scalar  function  introduced  in  Section  4.5.  The  Gaussian  function  has  been  chosen 
to  check  the  output  of  the  ADS  routine.  The  initial  domain  is  square:  [—0.5, 1.5]  x  [—0.5, 1.5] 
and  the  parameters  of  the  Gaussian  are  /x  =  (0.5  0.5)  and  a1  =  diag(0.1,  0.01).  The  expansion 
order  is  chosen  to  be  10,  the  error  tolerance  10-5  and  the  maximum  number  of  splits  10.  As  can 
be  seen  in  Figure  5.10  the  ADS  routine  applied  to  the  Gauss  function  returns  a  manifold  of  64 
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X 


X 


Figure  5.8:  Exact  solution  vs.  ADS  Manifold  results  for  sine  function  with  an  expansion  order  5 
and  tolerance  10-4.  On  the  left  side  there  is  the  exact  solution  on  the  single  domain  expansion, 
on  the  right  side  the  manifold  representation.  The  black  dots  on  the  x-axis  define  the  extremes 
of  the  sub-domains. 


Figure  5.9:  Representation  of  the  inverse-square  function.  On  the  left,  the  exact  solution  and 
the  single  domain  expansion,  on  the  right,  the  manifold  representation  and  the  extremes  of  the 
obtained  sub-domains 

Patches.  To  appreciate  the  accuracy  of  the  expansion,  the  difference  between  the  actual  values 
and  the  polynomial  evaluation  is  shown  in  Figure  5.11  for  different  maximum  number  of  splits. 
Figure  5.11(a)  is  the  result  of  the  evaluation  with  a  single  domain.  Figure  5.11(b)  was  obtained 
with  a  limit  of  4  splits  for  each  sub-domain  and  Figure  5.11(c)  is  the  result  for  the  limit  of  10 
splits  for  each  sub-domain.  To  understand  the  link  between  the  number  of  sub-domains,  the 
variation  of  expansion  order  and  the  tolerance.  Figure  5.12  is  shown.  The  results  are  acheived  by 
exploiting  the  Gaussian  function.  The  application  of  the  ADS  routine  to  the  Gaussian  function 
can  be  useful  the  describe  the  Patch  structure  and  its  similarity  with  the  chart  of  the  manifold. 
In  Figure  5.13  one  element  of  the  Patch  is  considered  and  is  mapped  into  the  DA  Manifold. 
The  Patch  has  a  chart  structure  so  it  must  contain  the  domain  and  its  map  onto  the  manifold. 
The  domain  is  defined  by  the  SplittingHistory  which  stores  the  story  of  the  domain  splits  and 
function  that  maps  the  sub-domain  onto  the  manifold. 

In  this  way  it  is  possible  to  have  a  complete  and  lightweight  structure  of  the  manifold  with  all 
the  information  about  the  resulting  domains  and  maps  after  the  ADS  routine. 
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Figure  5.10:  Gauss  function  exact  solution(left)  vs.  manifold  representation(right) 


Figure  5.11:  The  error  between  the  real  values  and  the  polynomial  approximation  is  computed 
for  the  single  domain  expansion  (a),  the  manifold  of  14  sub-domains  (b)  and  a  manifold  with  64 
sub-domains  (c) 


Figure  5.12:  Link  between  number  of  subdomains,  order  and  tolerance  for  Gaussian  function 

5.2.3  Multi  Variable  Vector  Function 

Recalling  the  definition  in  Section  4.6,  these  are  functions  of  the  type  /  :  Recalling 

the  goal  of  this  work,  two  main  test  cases  are  here  analyzed:  the  conversion  of  orbital  state 
representation  and  the  DAIOD  function. 
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Patch 


Figure  5.13:  Structure  of  the  Patch  and  its  components 


5.2.3.1  Conversion  of  orbital  state  representation 


In  Astrodynamics,  sometime  one  state  representation  better  suits  some  situations  than  another 
and  thus  a  conversion  function  is  needed.  It  is  useful  to  create  these  functions  in  the  DA 
environment  in  order  to  use  them  in  a  DA  domain  representation  to  obtain  the  same  DA  domain 
into  another  representation.  The  algorithms  for  these  transformations  contain  mathematical 
operations  and  non-linear  transformations.  The  highly  non-linear  transformations  which  do 
not  converge  well  when  expanded  around  a  single  point  could  generate  a  resulting  excessive 
error  of  the  Taylor  expansion.  If  the  conversion  function  is  used  as  map  function  in  the  ADS 
routine,  one  can  convert  the  domain  and  set  a  threshold  for  the  maximum  error.  In  this  work 
the  following  conversions  are  implemented: 
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where  ( a,  e,  i,  n,  w,v)  are  the  KEP,  (  v>  f>  g>  K  f,  l)  are  the  Modified  Equinoctial  Element  (MEE) 
and  ( y>  v vy>  vz )  are  the  Cartesian  coordinates  in  ECI  reference  frame.  Details  about  the  state 
representations  and  the  conversions  are  discussed  in  more  detail  by  Battin  (1999),  Curtis  (2015) 
and  Bate  et  al.  (1971).  They  have  been  adapted  to  the  DA  environment  such  that  they  map  from 

3ft6  to  3ft6. 

Furthermore  the  different  state  representations  could  have  one  or  more  singularity  points  where 
the  elements  are  not  defined.  Each  of  state  representations  can  be  expressed  as  a  vector  of  6  DA 
variable  and  mapped  onto  the  desired  state  representation  to  obtain  the  associated  vector  of  6 
DA  variables,  each  one  containing  the  n-th  order  Taylor  expansion  of  the  conversion.  To  check 
the  conversion  error  for  singularity  points,  the  ADS  tools  analyses  the  convergence  of  the  maps 
and  manages  the  error  during  the  conversion. 

As  an  example,  a  LEO  object  is  considered  and  its  orbital  elements  defined.  The  ECI  was 
found  with  a  double  precision  transformation  of  the  KEPs.  As  shown  in  Table  5.1,  the  object 
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Table  5.1:  Definition  of  the  phase  space  vector  and  orbital  elements 


ECI 

]i  a 

KEP 

V 

X 

6578  km 

a 

8038.9  km 

y 

0  km 

e 

0.2 

z 

0  km 

i 

0.1  deg 

Vx 

-0.719  km/s 

n 

0  deg 

vy 

8.432  km/s  0.2 

CO 

30  deg 

Vz 

0.015  km/s  0.2 

V 

330  deg 

considered  lies  close  to  the  equatorial  plane  and  its  ascending  node  belongs  to  the  x-axis  of  ECI 
reference  frame.  Without  loss  of  generality,  only  vy  and  vz  are  considered  as  DA  variables.After 
choosing  the  2-nd  order  expansion,  10-5  order  threshold  and  maximum  15  splits,  the  resulting 
sub-domains  obtained  during  the  conversion  are  shown  in  Figure  5.14. 

The  fundamental  result  is  that  the  ADS  executes  a  big  number  of  splits  close  to  vz  =  0,  thus 
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Figure  5.14:  Resulting  phase  space  domains  on  the  vy-vz  plane  for  the  conversion  from  ECI  to 
KEP,  the  red  box  is  the  nominal  value 

when  the  orbit  becomes  equatorial  and  some  keplerian  elements  are  not  defined.  Due  to  the 
singularity  is  not  possible  to  achieve  the  conversion  for  the  value  vz  =  0  but  the  ADS  reaches  a 
good  approximation  by  increasing  the  number  of  Patches  close  to  the  undetermined  points.  To 
deduce  supplemental  results,  let's  consider  a  single  perturbation  of  the  state  vector  along  the  vz 
direction  and  a  2-nd  order  expansion.  The  trivial  variation  of  vz  could  cause  a  modification  of 
the  eccentricity  of  the  resulting  orbit  but  as  it  can  be  seen  in  the  left  side  of  Figure  5.16  and  5.18 
that  variation  is  small. 
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Figure  5.15:  Right  ascension(upper)  and  inclination(bottom)  with  respect  to  the  perturbation  of 
vz  around  its  mean  value(red  square)  and  error  threshold  10-2 


Vz[km/s]  Vz[km/s] 

Figure  5.16:  Kepler  element  eccentricity(left)  and  absolute  error(right)  with  respect  to  the 
perturbation  of  vz  around  its  mean  value(red  square)  and  error  threshold  10-2 


In  Figure  5.15  and  5.17,  regardless  of  the  error  tolerance  chosen  for  the  ADS,  the  presence  of  a 
singular  point  in  the  neighborhood  of  the  mean  value  leads  to  the  indetermination  of  the  right 
ascension.  By  reducing  the  error  tolerance  the  Taylor  expansion  precision  increases  as  shown 
on  the  right  side  of  Figure  5.16  and  5.18.  The  absolute  error  is  computed  at  the  extremes  of  the 
domain  space  between  the  actual  values  and  the  evaluation  of  the  KEP  map.  The  reduction  of  the 
tolerance  produces  a  decrease  in  error  magnitude.  Regardless  of  the  reduction,  the  singular  point 
corresponding  to  vz  —  0  does  not  allow  the  Taylor  expansion  to  approximate  well  enough  the 
Keplerian  element,  thus  next  to  this  point  there  is  a  peak  of  the  absolute  error  which  produces  a 
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consecutive  split  of  the  domains.  The  domains  of  minimum  size  are  located  in  the  neighborhood 
of  the  singular  point  as  shown  on  the  left  side  of  Figure  5.16  and  5.18. 


Vz[km/s\ 


Figure  5.17:  Right  ascension(upper)  and  inclination(bottom)  with  respect  to  the  perturbation  of 
vz  around  its  mean  value(red  square)  and  error  threshold  10-4 


Figure  5.18:  Kepler  element  eccentricity(left)  and  absolute  error(right)  with  respect  to  the 
perturbation  of  vz  around  its  mean  value(red  square)  and  error  threshold  10-4 


41 

DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


5.2.3.2  Initial  orbit  determination  (IOD)  and  Automatic  domain  splitting  (ADS) 

In  Section  5.2.3  another  multi-variable  vector  function  of  interest  was  introduced,  which  corre¬ 
spond  to  the  solution  of  the  Initial  orbit  determination  (IOD).  In  this  situation,  the  map  is  the 
IOD  algorithm  itself,  defined  in  Section  3.3.  The  IOD  function  takes  as  input  the  Observation 
Set  and  returns  the  TPS  of  the  object  state  at  the  central  time  of  the  observation  i.e  (*2,^2)- 
The  input  is  detailed  in  Section  1.1.  The  map  is  a  function  /  :  x  e  ►  f{x)  e  3?^,  with  the 
subscript  being  the  expansion  order  of  the  TPS  and  it  can  be  analysed  by  the  ADS  to  control 
the  error.  The  ADS  needs  to  be  slightly  adjusted  to  be  able  to  get  as  input  the  IOD  algorithm 
and  all  the  parameters  needed  for  it:  indeed,  a  structure  parameter  is  also  passed  to  the  routine. 
It  contains  the  angular  observations  [ol\,  5\,  0L2,  $2,  ^3,  ^3),  the  three  observation  times,  the  three 
observation  precisions  (cp/L,  crP/2/  (7p/3)  and  the  observatory  position  at  the  three  observation 
times  obtained  by  using  SPICE  as  described  in  Section  2.1. 

The  ADS  routine  with  the  multi-variable  vector  function  can  be  validated  by  choosing  the  orbit 
of  a  real  object  and  using  the  VO  to  generate  the  observation.  Given  the  coordinates  of  a  satellite 
in  term  of  its  TLE  (see  Table  5.2),  six  simulated  Observation  Sets  were  generated  with  the  VO, 
each  with  a  different  time  interval  between  the  single  observations. 

As  a  first  step,  a  6-th  order  expansions  was  performed  and  the  sensitivity  of  the  results  to  the 

Table  5.2:  Satellite  TLE 

0  CADRE 

1  41475U  98067HV  16172.13123003  .00037064  00000-0  49179-3  0  9994 

2  41475  51.6443  42.8436  0003066  161.2134  270.8982  15.57867415  5398 

error  tolerance  was  analysed.  In  Figure  5.19  the  absolute  errors  for  the  different  Observation 
Sets  are  shown.  It  is  evident  that  the  set  with  a  longer  time  interval  between  two  consecutive 
observations  generate  a  TPS  object  state  with  higher  precision.  In  Figure  5.20  the  number  of 
sub-domains  generated  by  the  ADS  is  shown  and  it  can  be  seen  that  the  number  of  sub-domains 
starts  increasing  at  relatively  low  tolerances,  i.e.  when  the  absolute  error  starts  decreasing  (see 
Figure  5.19). 

In  a  second  analysis,  the  error  tolerance  is  fixed  to  10-5  and  the  expansion  order  varied.  Figure 
5.21  shows  the  variation  of  the  absolute  error.  It  is  clear  that  the  error  decreases  by  raising  the 
order  but  it  is  also  possible  to  see  in  Figure  5.22  that  the  lowest  expansion  order  is  not  accurate 
enough  to  be  able  to  estimate  the  error,  thus  the  ADS  does  not  split  the  initial  domain,  even 
though  the  absolute  error  is  beyond  the  tolerance. 


5.2.3.3  Admissible  region  (AR)  and  automatic  domain  splitting  (ADS) 

The  ADS  can  also  be  used  to  increase  the  precision  of  the  AR  obtained  in  Section  5.1.  Now  the 
first  two  ARs  obtained  are  re-evaluated  with  the  ADS  routine:  the  accurate  GEO  observation 
and  the  inaccurate  GEO  observation.  The  ADS  routine  is  started  with  a  3-rd  order  expansion 
for  the  TPS.  Figures  5.23  and  5.24  refer  to  the  accurate  observation  while  Figures  5.25  and 
5.26  refer  to  the  inaccurate  observation.  The  black  box  on  top  of  Figures  5.23  and  5.25  is  the 
box-shaped  AR  obtained  in  Section  5.1,  while  the  boxes  shown  in  Figures  5.24  and  5.26  are  the 
ARs  found  for  each  of  the  sub-domains  created  by  the  ADS  routine.  The  bottom  plot  of  each 
figure  represents  the  evaluation  of  the  ARs:  the  red  dot  is  the  TPS  evaluation  while  the  blue 
triangles  are  the  pointwise  values  for  the  range  and  range-rate,  which  are  computed  by  solving 
the  IOD  algorithm  in  the  variations  5oc  and  55.  Comparing  these  values  it  is  possible  to  see  that 
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Figure  5.19:  Absolute  Error  for  the  set  of  observations  with  different  time  interval,  6-th  order 
expansion  and  varying  error  tolerance 
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Figure  5.20:  Number  of  sub-domains  for  the  set  of  observations  with  different  time  interval,  6-th 
order  expansion  and  varying  error  tolerance 
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Figure  5.21:  Results  of  Absolute  Error  for  the  set  of  observation  with  different  time  interval,  10  5 
error  tolerance  and  varying  expansion  order 
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Figure  5.22:  Number  of  sub-domains  for  the  set  of  observations  with  different  time  interval, 
10-5  error  tolerance  and  varying  expansion  order 
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the  differences  between  the  actual  value  and  the  single  domain  evaluation  one  are  large,  indeed 
the  ADS  splits  the  initial  domain  and  it  generates  a  manifold,  which  then  well  describes  the  AR, 
given  the  similarity  of  the  results.  Furthermore  this  result  allows  to  further  reduce  the  size  the 
first  AR  obtained  in  Section  5.1,  which  was  already  significantly  smaller  than  the  AR  computed 
from  the  literature. 

At  the  end  the  analysis  of  the  IOD  performed  with  the  ADS  it  can  be  noted  that  the  sub- 
domain  bounds  for  the  accurate  observation  are  uniform  within  the  AR  obtained.  This  result 
is  connected  to  the  orbital  elements  of  the  OS  space  shown  in  Figure  5.27:  since  they  are  far 
from  the  singular  point,  the  ADS  routine  works  to  sample  the  AR  and  then  checks  the  error 
to  achieve  higher  precision  with  the  IOD  algorithm  quite  smoothly.  On  the  other  hand,  the 
inaccurate  observation  places  the  point  solution  close  to  the  singular  elements,  as  it  can  be 
seen  for  the  eccentricity  in  Figure  5.28.  Hence  the  ADS  routine  needs  to  check  the  error  in 
the  TPS  but  also  has  to  manage  the  impossibility  to  determine  the  OS  in  the  singular  point. 
To  support  this  assertion,  it  is  possible  to  see  the  large  number  of  sub-domains  close  to  the 
value  of  p  that  corresponds  to  the  singularity  on  e  and  co  (compared  Figure  5.26  with  Figure  5.28). 
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Figure  5.23:  AR  and  ADS  results,  single  domain  evaluation  vs.  actual  values.  The  big  black 
box-shape  is  the  AR,  the  black  square  is  the  mean  value  of  the  AR,  the  red  dots  are  the  single 
domain  evaluations  and  the  blue  triangles  are  the  actual  object  state. 

- Initial  bound  evaluation  *  Initial  domain  actual  value  •  After-Split  manifold  evaluation  □  Nominal  Value 
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Figure  5.24:  AR  and  ADS  results,  manifold  evaluation  vs.  actual  values.  The  black  box-shape 
are  the  bounds  of  subdomains  obtained  by  ADS  are  add.  Accurate  observation. 
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- Initial  bound  evaluation  *  Initial  domain  actual  value  •  Initial  domain  evaluation  □  Nominal  Value 
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Figure  5.25:  AR  and  ADS  results,  single  domain  evaluation  vs.  actual  values.  The  big  black 
box-shape  is  the  AR,  the  black  square  is  the  mean  value  of  the  AR,  the  red  dots  are  the  single 
domain  evaluations  and  the  blue  triangles  are  the  actual  object  state. 
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Figure  5.26:  AR  and  ADS  results,  manifold  evaluation  vs.  actual  values.  The  black  box-shape 
are  the  bounds  of  subdomains  obtained  by  ADS  are  add.  Inaccurate  observation. 
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Figure  5.27:  Orbital  element  of  the  corresponding  actual  object  state  of  Figure  5.25.  Accurate 
observation. 


Figure  5.28:  Orbital  element  of  the  corresponding  actual  object  state  of  Figure  5.23.  Inaccurate 
observation 
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Chapter  6 

Conclusion 


This  report  has  detailed  the  work  of  the  team  in  the  period  June,  16  2016  -  January,  16  2017.  The 
Universidad  de  la  Rioja  (UR)  team  has  completed  the  development  of  the  IOD  algorithm  and 
it  was  implemented  for  the  optical  observation.  This  algorithm  was  developed  by  exploiting 
the  well-known  Gauss's  and  Lambert's  problem  and  translating  them  in  theDA  environment. 
To  validate  the  IOD  algorithm  the  VO  was  used.  The  VO  gets  the  TLEs  of  real  known  objects 
and  generates  the  Observation  Set  for  a  given  At  as  detailed  in  Section  2.2.  The  results  obtained 
by  the  IOD  were  analyzed  in  terms  of  the  AR,  already  known  in  literature.  The  model  was 
compared  with  the  previous  approaches  as  detailed  in  Section  5.1,  and  it  showed  a  significant 
improvement  in  the  description  of  the  uncertainty  area  associated  with  the  IOD  solution. 

The  PoliMi  team  focused  on  the  DA  representation  of  the  domain,  the  TPS  evaluation  and  the 
domain  propagation  in  the  DA  environment.  In  particular  the  DA  representation  of  a  mathe¬ 
matical  manifold  has  been  implemented,  with  the  creation  of  the  DAmanifold  as  highlighted  in 
Section  4.1  and  4.2.  The  DA  manifold  is  useful  to  represent  a  domain  and  its  map.  After  defining 
the  DA  Manifold,  it  was  possible  to  implement  the  ADS  tool  with  the  new  DA  structure,  thus 
obtaining  a  better  description  of  the  linkage  between  the  domains  and  their  maps  during  the 
sampling.  The  new  structure  was  implemented  for  the  single  and  multi-variable  scalar  function 
and  for  the  multi-variable  vector  function.  Thanks  to  the  improvements  introduced,  a  single 
routine  is  able  to  manage  these  three  different  cases.  Furthermore,  a  set  of  functions  to  transform 
between  different  orbital  state  representation  was  implemented  in  the  DA  framework. 

The  last  part  of  the  activity  focused  on  merging  the  IOD  and  the  ADS  routines.  The  objective  was 
to  obtain  an  algorithm  that  allows  for  suitably  controlling  the  error  in  the  IOD  thanks  to  the  ADS 
tool.  The  results  with  the  new  merging  model  were  compared  with  the  AR  from  the  literature 
and  with  the  AR  on  a  single  domain.  As  highlighted  in  Section  5.23.2,  the  IOD  performed  with 
the  ADS  routine  is  able  to  obtained  a  reduction  of  the  AR,  thus  a  more  accurate  OS,  both  for 
accurate  and  inaccurate  observations. 
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List  of  Acronyms 


ADS 

automatic  domain  splitting 

AIDA 

Accurate  Integrator  for  Debris  Analysis 

AOD 

accurate  orbit  determination 

AR 

admissible  region 

CAR 

constrained  admissible  region 

DA 

differential  algebra 

DAIOD 

Differential  Algebraic  Initial  Orbit  Determination 

DAvector 

differential  algebraic  vector 

ECI 

Earth-centered  inertial 

FOV 

field  of  view 

GEO 

geostationary  Earth  orbit 

GTO 

geostationary  transfer  orbit 

HEO 

high-Earth  orbit 

IOD 

initial  orbit  determination 

KEP 

Keplerian  classical  Element 

LEO 

low-Earth  orbit 

MEE 

Modified  Equinoctial  Element 

MEO 

medium-Earth  orbit 

NAIF 

Navigation  and  Ancillary  Information  Facility 

NASA 

National  Aeronautics  and  Space  Administration 

OD 

orbit  determination 

ODE 

ordinary  differential  equation 

OS 

orbit  set 

PDF 

probability  density  function 

SA 

short  arc 

SPICE 

Spacecraft  Planet  Instrument  C-matrix  Events 

TLE 

two-line  elements 

TPS 

truncated  power  series 

52 

DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


TSA 

UR 

VO 

VSA 


too  short  arc 
Universidad  de  la  Rioja 
virtual  observatory 
very  short  arc 
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